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(*~^ ' We present molecular dynamics simulations of mono- or bidisperse inelastic granular gases driven 

' by vibrating walls, in two dimensions (without gravity). Because of the energy injection at the 

^SJ . boundaries, a situation often met experimentally, density and temperature fields display heteroge- 

neous profiles in the direction perpendicular to the walls. A general equation of state for an arbitrary 
mixture of fluidized inelastic hard spheres is derived and successfully tested against numerical data. 
Single-particle velocity distribution functions with non-Gaussian features are also obtained, and the 
influence of various parameters (inelasticity coefBcients, density. . . ) analyzed. The validity of a 
recently proposed Random Restitution Coefficient model is assessed through the study of projected 
collisions onto the direction perpendicular to that of energy injection. For the binary mixture, the 
non-equipartition of translational kinetic energy is studied and compared both to experimental data 
i-C ' and to the case of homogeneous energy injection ("stochastic thermostat"). The rescaled velocity 

^ ' distribution functions are found to be very similar for both species. 

e: 

I. INTRODUCTION 
> , 

, Due to the intrinsic dissipative character of inter-particle collisions, an energy supply is requested to fluidize a 

■ granular gas. This is often achieved by a vibrating boundary, and the resulting vibro-fluidized beds provide non trivial 
realizations of non equilibrium steady states. The understanding of such far from equilibrium systems requires a correct 
description of the energy exchange between the vibrating piston and the granular medium, as well as a macroscopic 
continuum theory to describe the evolution of the relevant coarse-grained fields ^ (density, temperature etc.). In 

(~| particular, the derivation of an accurate equation of state is a key step in the hydrodynamic approach. 

O A simple, fair and much studied theoretical framework to capture the inelastic nature of grain-grain collisions in 

, O, . a rapid granular flow is provided by the inelastic hard sphere model In this article, we present the results 

of extensive molecular dynamics (MD) simulations of inelastic hard spheres driven by an energy injection at the 
T-H , boundaries, for both a one component fluid (mono-disperse case) and a binary mixture (bidisperse situation). We 

^ ' analyze in detail the effects of several parameters that may be difficult to tune experimentally, with a particular 
C [ emphasis on the profiles of the hydrodynamic fields. 

■ This article is organized as follows: in section II, we present the model and derive an equation of state for an 
' arbitrary mixture of inelastic hard spheres, going beyond the ideal gas contribution in view of performing accurate 

hydrodynamic tests. The equation of state obtained is a natural generalization of its standard counterpart for elastic 
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hard spheres. The two following sections ([II and IV) are then devoted to molecular dynamics simulations for one- 



component systems and for binary mixtures. In both cases, we restrict ourselves to two-dimensional simulations, both 
for simplicity and for comparisons with 2D experimental data ^, 0, As in the experiments, the energy loss due to 
[ inelastic collisions is compensated for by an energy injection by vibrating or thermal walls, which leads to heterogeneous 
density and temperature profiles. The various profiles and velocity distribution functions are studied and compared 
with experiments whenever possible. Moreover, projecting the dynamics onto the direction perpendicular to that of 
energy injection allows to assess the validity of the random restitution coefficient model proposed in [lO|. The 



influence of various parameters on the non-equipartition of energy in a binary mixture is studied in section [V , and 
^ , comparison with experimental data and with the case of homogeneous energy injection is performed. In this latter 
case, the velocity distribution functions are analyzed and shown to be very similar for the two species. Conclusions 
are finally presented in section 



II. THE MODEL - COMPUTATION OF AN EQUATION OF STATE 

We consider a mixture of Afs species of hard spheres in dimension d, with diameters ai and masses rrii, where 
1 < i < Mg. A binary collision between grains of species i and j is momentum conserving and dissipates kinetic 
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energy. In the simplest version of the model, the collision i-j is characterized by one inelasticity parameter: the 
coefficient of normal restitution aij . Accordingly, the pre-coUisional velocities (vi , Vj ) are transformed into the post- 
collisional couple {vl,Vj) such that 

- aij){a ■ v.,j)a (1) 



TO," + Too 



rui + rrij 

where Vij — — Vj and a is the center to center unit vector from particle i to j. Note that — cxji to ensure the 
conservation of total linear momentum rriiVi + nijVj. 

We also considered an extension of the previous model allowing for rotations, introducing a coefficient of tangential 
restitution a*^- Q, see appendix]^. The coUision law (0)-(||) is then recovered for a-^ = — 1. 

Irrespective of the value of the tangential restitution coefficient a* , the linear- momentum change for particle i in a 
collision i-j reads 

5pi = -6pj ^ {l + a^j){a ■Vij)a. (3) 

mi + rrij 

In appendix |b[ we use this relation to compute an equation of state for the homogeneous isotropic mixture, invoking 
the virial theorem (the pressure is defined kinetically from the momentum transfer and does not follow from a 
statistical mechanics derivation). The total density is denoted p and the partial densities pi — Xip (the number 
fractions Xi are such that '^^Xi — 1). The temperature of species i is Ti, defined from the mean kinetic energy 
of subpopulation i: dTi — {rriiV^. Only for an elastic system is the energy equipartition Ti = T,\fi recovered 
H, |l2|, |l3|, 1^, 1^, 111. It is found in appendix g that the pressure in dimension d reads 

P = J2p^T, + prj2'-' Y.^,x, ^^^^ (1 + a.,) T, ^ (4) 

independently of q;*^, where (Ty = [pi + o'j)/2, (cr'^) = J^i^i'^t^ V *he packing fraction (e.g. t] = 7rp(cr^)/6 in 
three dimensions), and the Xij ^re the pair correlation functions at contact. The latter -unknown- quantities may be 
approximated by their elastic counterparts (see for a general procedure to infer reliable pair correlation functions 
in a multi-component d-dimcnsional hard-sphere fluid from the equation of state of the mono-disperse system). In 
the following analysis, it will turn sufficient to include only the low density behaviour Xij = 1 to improve upon the 
ideal equation of state P = P^'^^^^ = J2iPi'^i> ^^^^ holds for p ^ only. We emphasize that no approximation has 
been made on the single-particle velocity distribution in the derivation of Eq. (Q) (the key assumption is that the 
two-body distribution function factorizes at contact in a product of the single-particle distribution ) . 

It is instructive to check the validity of our equation of state by considering the elastic limit where aij = 1 and 
Ti — T. A straightforward calculation (under the reasonable and often made assumption that Xij — Xji) shows that 
the mass ratio simplifies and expression (^) may be cast in the form 

^ = l-Hr;2'^-i^2;,x, -^X.„ (5) 

which is the correct result (see e.g. [251 ). In particular, for the single species (mono-disperse) problem, we recover 
the virial relation P/{pT) = 1-1- 2'^~^rix- 

We finally generalize Eq. (^) to the situation of a continuous size distribution, with a probability density distribution 
W{(j) (normalized to 1 so that ((t") — J a^'^W); the temperature is in general a continuous function T{a) of size and 

- = / Wia)Tia)da + ^ f da da' W {a)W {a' ) (1 + a..O T{a) ^^^^^ X..' ■ (6) 

P J 2 J TO^ m^' (cr") 

In the following sections, the above equation of state will be used to test hydrodynamic predictions for a monodisperse 
system and for a binary mixture. 



III. MOLECULAR DYNAMICS SIMULATIONS FOR THE ONE COMPONENT SYSTEM 



A. Introduction 



We have implemented molecular dynamics simulations with an event-driven algorithm for N spherical particles in a 
two-dimensional i x i box. Periodic boundary conditions are enforced in the x direction, while the energy loss due to 
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collisions is compensated by an energy injection by two walls situated at y = and y = L (we consider the amplitude 
of motion of the walls to be small so that their positions are considered as fixed , which avoids the complication of 
heat pulses propagating through the system |2^). We will refer to the y direction as the "vertical" one, although we 
are interested in regimes for which gravity can be neglected |^ (i.e. when the shaking is violent enough). The energy 
can be injected in two ways: 

• by thermal walls which impose a given temperature of order Tq |p7| : when a particle collides with the 
wall, its new vertical (along y) velocity is extracted at random from the probability distribution function 
w/\/7^exp[— w^/(2-yTo)], whereas Vx is unaffected. 

• by vibrating walls: for simplicity, we consider walls of infinite mass moving in a sawtooth manner: all particles 
colliding with a wall find it with the same velocity > at y = 0, —vq at y = L. The particle-wall collisions 
are considered elastic. A particle of velocity v with Vy < colliding with the bottom wall at i/ = (resp. Vy > 
at the upper wall) sees its velocity change to v' according to Vy = 2vo — Vy (resp. v'y ~ —2vq — Vy), whereas the 
x-component is unaffected {v'^ = Vx)- 

In both cases, energy is injected in the vertical direction only, and transferred to the other degrees of freedom through 
inter-particle collisions. The vibrating walls being the situation closer to the experimental one, most of our results 
will be presented in this case, and the effect of injection modes will be briefly discussed. 

In this section, we consider the monodisperse case: all particles have the same mass m[— 1), diameter cr, restitution 
coefficients a and a*. Most of the simulations are done with N = 500 particles, some with N = 1000 particles (low 
enough to avoid clustering or inelastic collapse) . For our two-dimensional system, the local packing fraction at height 
y, where the local density is p{y), is defined as ri{y) = tt p{y)a'^ / A. The global (mean) packing fraction is denoted 770: 

^0 = JoViy)dy/L. 

Starting from a random configuration of the particles (with the constraint of no overlap), we let the system evolve 
until a steady state is reached. Data on density and temperature profiles as well as on velocity distributions are 
monitored as a time averages; the various quantities are averaged along the x direction since the system remains 
homogeneous in this direction. 



B. Density and temperature profiles 

The first observations concern the density and temperature profiles: Figs. |] and ^ show that the density is lower 
near the walls, where the temperature is higher as expected since energy is injected at the walls and dissipated in 
the bulk of the system The profiles are qualitatively similar for thermal or vibrating boundaries. Moreover, the 
whole temperature profile is proportional to the temperature Tq imposed by a thermal wall or to the square of the 
velocity vq of the vibrating boundary, while a change in Tq or vq does not change the density profile (not shown). 
As the mean density increases or a decreases, the profiles get more heterogeneous; as a* is increased, more energy 
is transferred to rotational degrees of freedom, so that the temperature decreases, while the density profiles become 
slightly more peaked (Fig. ||). 

Fig. ^ clearly shows another feature resulting from the energy injection into the vertical direction: the temperature 
is anisotropic, i.e. (v^) ^ (Vy), with Ty > T > Tx- The anisotropy A{y) = {Ty — Tx)/{2T) is larger at the boundaries, 
where energy is fed into the vertical direction, decreases due to inter-particle collisions, and reaches a plateau in the 
middle of the slab. The plateau value decreases for increasing number of particles or increasing densities (not shown), 
as in experiments m ; the global anisotropy profile and the plateau values are comparable to experimental values m . 



C. Equation of state and hydrodynamics 

The equation of state derived in section II reduces, in the case of a two dimensional one-component homogeneous 
system, to the relation 

pT[l + il + a)r,x], (7) 

where Xj the pair correlation function at contact, depends on the packing fraction r). We will use the form x = 
(1 — 777/16)7(1 — 77)^, which has been shown to be accurate for elastic hard disk liquids p9t . 

The Irydrodynamic equations (see appendix |C|, and ||l|) lead to dyP = in the absence of a flow field. We check 
in Fig. ^(a) the constancy of P with y by plotting the ideal gas contribution p{y)T{y) (lines) and P{y) given by Eq. 
(0) (i.e. ideal gas contribution plus Enskog correction). While, at small enough densities (not shown), p{y)T{y) is 
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constant in the bulk (y G [0.2L, 0.8L]), the Enskog correction is necessary for the densities used in Fig. ^ (note that 
the density can be quite larger in the middle of the system than at the boundaries). We also note that the inelasticity 
term (1 + a) is relevant [the profiles of pT (1 + 2rix), not shown, do not display a uniform shape with y]. In all cases, 
boundary layers (y < 0.2L and y > 0.8i) are observed in which the pressure decreases. This discrepancy can be 
related to the anisotropy described in the previous subsection (pressure and temperature are most anisotropic near 
the walls). 

The comparison with hydrodynamics may be improved as follows. The pressure tensor P is diagonal in the present 
no flow situation, but has different xx and yy components, and the homogeneity along the x direction implies that the 
condition of vanishing flow field V P — reduces to dyPyy = 0. We therefore check in Fig. ^b) that the yy component 
of the pressure tensor, given by the equation of state ^ with the total temperature T — [T^ + Ty)/2 replaced by 
its vertical component T^, is uniform in the whole system. With Enskog correction, the corresponding profiles are 
remarkably flat. This result could be tested in experimental situations in which both and Ty are measured. Such 
an analysis validates both the hydrodynamic picture and the equation of state proposed by automatically sampling 
several densities in a single run. 

At low densities, assuming the ideal gas equation of state to hold, the hydrodynamic study of Rcf. (recalled in 
appendix]^), leads to the following analytical prediction for the temperature profile: 

y_ 

L 

e 

where Tq is the temperature at the boundaries and S^„i is proportional to the total number of particles. The cor- 
responding fits of the temperature profiles are shown in Fig. |5|; a good agreement is obtained, especially at lower 
densities as expected [since the ideal gas equation of state is a crucial ingredient in the derivation of Eqs. (H)]. We 
use one fitting parameter to obtain T/Tq. Fig. ^showed that consideration of the "vertical" pressure Pyy led to a 
better agreement with hydrodynamic predictions than the total P^x + Pyy A similar conclusion is incorrect for the 
temperature profiles: the transport equation for the temperature is scalar [see Eq. 
total T, not for the vertical Ty. 

D. Velocity distributions 

Because of the energy injection through the walls, the velocity distributions are anisotropic, and a priori depend on 
the distance to the walls. The vertical velocity distribution also depends on the nature of the walls as shown in Fig. ||. 
A smooth distribution is obtained in the vicinity of a thermal wall, while the incoming and out-coming particles yield 
two separated peaks for vibrating walls (see also [|l|). 

On the other hand, the reseated horizontal velocity distribution P{cx) (with — Vx/^/T^) is remarkably independent 
of the distance from the walls (outside the boundary layers), even if the temperature changes with y. Fig. |^ shows 
clearly non-Gaussian features similar to the experimentally observed ones [^, ^ |3^ , with in particular overpopulated 
both small-velocity and high-velocity regions. A slight dependence on the parameters is obtained: P{cx) broadens if 
the inelasticity increases (i.e. if a decreases), if a* increases, or if 770 or N increase. Experimentally, the dependence 
on density or material properties is weak and difficult to measure |q| but seems to exist, in particular as far as N is 
varied The angular velocity distributions, also displayed in Fig.^ share a similar non-Gaussian character and the 
same dependence with the parameters. 

As density or inelasticity are further increased, clustering phenomena may occur, leading to heterogeneities along 
the X direction, with the coexistence of colder, denser regions with hotter, less dense ones. The average over the x 
direction then leads to artificially broad P{cx). 

Finally, as a general rule, thermal walls lead to slightly broader velocity distributions than vibrating walls. 

E. Effective restitution coefficients 

We now turn to the study of the effective inelasticities introduced in the context of a Random Restitution Coefficient 
model (RRC) [|, |l§: even if the restitution coefhcient a is constant, the energy is injected in the vertical direction 
and transferred to other degrees of freedom through collisions, so that the ejjeetive restitution coefficient for collisions 



g + sinhgcosh(^,n - £,) 
^rn + sinh^m 




(8) 



C2)] and Eqs. (|) hold for the 
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projected onto the x direction, 



aid = — , (9) 

Vl2,x 



may be either smaller or larger than 1. This leads to the definition of a one-dimensional effective model with a 
restitution coefhcient taken at random from a given distribution at each collision p^ . 

Values of aid have been experimentally measured |l^ and shown to display a broad probability distribution 
fi{aid) very similar for various materials and densities. We have measured aid for many collisions and thus obtained 
its distribution, displayed in Fig. |^ together with experimental data for steel and glass beads. A remarkable agreement 
is found. Our study shows that ii{aid) display a a^^ tail for aid > 1, irrespective of a, a* and density. 

The importance of the correlations between aid and the relative velocity g — Vyij^f^ of the colliding particles 
has been emphasized in [l0[| and is revealed by the computation of p^{aid\g^ , the distribution of aid conditioned by 
a given value of g^', although no precise experimental determination of the conditional \x{aid\gx) could be achieved 
in pQ l, strong evidences for a sharp cut-off oc ^jgx at large values of aid were provided and the form [i,(aid\gx) oc 
exp (^(aici(73;)^/i?) at large aid has been proposed. The conditional \i{aid\gx) obtained in the present MD simulations 
confirm the above picture; they are displayed in Fig. ^ and show an exp [— {aidgxf' I decrease for the case of vibrating 
walls (closer to the experimental situation), and a broader form ex:^\—{aidgx)li^\ for thermal walls. Moreover, 
although [i{aid) is not sensitive to the various parameters, the cut-off K increases [i.e. leads to broader [i(aid\g3^\ if 
a decreases, and if a* or 770 increase. 

These findings, together with the evolution of the velocity distributions P{cx) with the parameters, is in com- 
plete agreement with the one-dimensional effective RRC model put forward in ilQ], for which broader conditional 
distributions [i(aid\gx) are linked to broader P(cx) (at large Cx, compared to the Gaussian). 

Finally, the energy restitution coefficient 

(10) 

may also be viewed as a random variable that can take values larger than unity due to energy transfers between 
rotational and translational degrees of freedom 10 . Fig. displays the p.d.f. p(/?) obtained in the MD simulations 



for various values of a , together with the experimental data of [g| |10| for steel beads. yo(/3) becomes wider as a is 
increased, but the experimental distribution is broader, which may be traced back to the fact that in the experiments 
mentioned above, the beads can rotate in three dimensions whereas our simulations are limited to 2D rotations. 

IV. MOLECULAR DYNAMICS SIMULATIONS FOR THE BINARY MIXTURE 

In this section, we investigate the properties of vibrated binary mixtures; such systems have recently attracted 
much attention, both on the experimental H |l|, |l|l and theoretical side ||, [l5[ |l|, [13, |l|, |l|, |l], H |3^, H . In 



particular the breakdown of energy equipartition between the two constituents of the mixture has been thoroughly 
investigated. 

The main difference with previous studies consists here in the realistic character of both MD simulations (as opposed 
to Monte Carlo methods) and the energy injection mechanism at the boundaries; the set-up is the same as in the 
previous section, with however two types of particles, with masses mi, TO2, sizes tri, 02- The three normal restitution 
coefficients (corresponding to the three possible types of collisions) are an, an = a2i, (222 ■ In the context of a 
forcing mechanism through a random external force |2j, |3^, it has been shown that the influence of size ratio on 
the temperature ratio measuring the energy non-equipartition was rather weak ]l^ compared to that of inelasticity 
parameters or mass ratio. We shall consequently limit our study to identical sizes CTi = (T2 in two dimensions, which 
corresponds to the experimental situation we will refer to |^ ^. For simplicity, the tangential restitution coefficients 
alj are also taken equal. 

As in the monodisperse case, we measure density and temperature profiles, velocity distributions as well as the 
temperature ratios 7(y) = T2{y)/Ti{y), ^x{y) = T2,x{y) /Ti^xiy)-, lyiv) = T2,y{y)/Ti^y{y). Some comparison with 
experimental data g] will be proposed whenever possible. 

A. Equation of state 

We first test the equation of state (||) in Fig. |ll|. As in the monodisperse case, the Enskog correction is clearly 
relevant, even at low global densities, since the density profiles reach relatively high values for y ~ L/2. It is however 
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sufficient to truncate the equation of state at second virial order, which amounts to take the low density Hmiting value 
Xij = 1 for the pair correlation functions at contact: 

9 

P ~ piTi + P2T2 + — ■ [(1 + aii)plm2Ti + (1 + ai2)/5i/02(™iT'2 + "I2T1) + (1 + a22)p2'^iT'2] . (11) 

Z(mi + 1712) 

Moreover, the boundary layer in which the anisotropy is strong is still apparent if the global temperatures Ti and T2 
are used, while use of the vertical ones Ti^y and TJj^y, suggested by the anisotropy of temperatures and pressure as in 
the monodisperse case, leads to a uniform yy component of the pressure tensor in the whole system. The functional 
dependence of pressure upon densities is therefore accurately reproduced by the equation of state (|Tl|) . 

Although we have not extended the hydrodynamic approach of Brey et al. ||l| to binary mixtures (it would be 
possible making use of the Navier-Stokes like equations derived in jSlj l where only the overall temperature associated 
with both species serves as a hydrodynamic field, but where the transport coefficients explicitly depend on temperature 
ratio), we see in Fig. |l^ that the temperature profiles can be fitted, at low density, by the form (^). We emphasize 
that there is no fundamental reason for the agreement. The quality of the fit is much better for the less massive 
particles whose density is more homogeneous across the system (see next subsection). For simplicity, we have used 
the short hand notation aij = 0.7; 0.8; 0.9 for the situation where an = 0.7, ai2 = 0.8 and 022 = 0.9. 



B. Non equipartition of translational kinetic energy 



The density and temperature profiles are displayed for various values of the parameters in Figs. 13 and |1J. The 
more massive particles (labeled 1), which display a more heterogeneous profile and are denser in the middle of the cell, 
have typically larger kinetic energies than the lighter ones: generically 7 = T2/T1 is smaller than 1, as in homogeneous 
mixtures [|l5[ ^ . The study of the y-dependence of 7 shows that 7 increases from the boundaries to the center of the 
system, and is constant across a wide range of y even if Ti and T2 vary significantly. As also experimentally shown 
in 7 is very close to 1 if mi = TO2j even if the inelasticities of the particles are different. It decreases if the mass 
ratio increases (Figs. ^), but displays only a very weak (but strikingly similar to experimental data) sensitivity on 



the global density (Fig. rsj) as well as on the relative densities of heavy and light particles; moreover, 7 may increase 
or decrease as ryi,o/'72.o is increased (see Fig. [l^ ), depending on the relative inelasticities. 

The anisotropy in the temperatures yield an anisotropic 7; we obtain, as in experiments y x > 7 > 7j,, with 
also different shapes: decreases from the walls to the center while 7 and jy increase (Fig.^Tsj). All these results 
are in very good agreement with the existing experimental results for two-dimensional vibrated mixtures |^ . We 
summarize in tables | and || some of the effects reported here. 

When rotations are included (and thus a* > — 1), 7 decreases. Moreover, the ratio of rotational kinetic energies 'jr 
can then be measured. As shown in tables | and 7^ takes values of the same order as 7. This quantity may also 
be computed from experimental data, although measures of rotational velocities are a priori more difficult than that 
of translational ones. 

The measured values of 7 are of the same order as the experimental data. We do not however try to obtain a 
precise numerical agreement for the following reasons: (i) in the experiments of the beads can rotate in three 
dimensions, whereas the simulated spheres rotate in two dimensions only; since a* has a strong effect on 7, we suspect 
that this difference between experiments and simulations may affect 7; moreover, the experimental value of a* is not 
known, and the precise validity of the inelastic hard sphere model with a tangential restitution coefficient should be 
assessed; (ii) different energy injection mechanisms (thermal vs. vibrating walls, homogeneous driving vs. injection 
at the boundaries) lead to different values of 7; even if the energy injection by vibrating walls is reasonably realistic, 
such a sensitivity of 7 renders its precise numerical prediction elusive. 

Nonetheless, the qualitative very good agreement, even for subtle effects (see e.g. Fig. |l5|), between numerics 
and experiments, and the possibility to change the various parameters in the simulations, allow us to make some 
predictions on the effect of various parameters: for example, increasing the mass ratio should yield smaller values 



of 7 (Fig. nM. Moreover, Fig. 14 makes it clear that the value of 7, at given mass ratio, is smaller for inelasticities 



aij — 0.9; 0.8; 0.7 than with "reverse" inelasticities Uij = 0.7; 0.8; 0.9. This effect was already noted in |16 and has 
the following intuitive interpretation: when the more massive particles are more inelastic, they loose more energy, 
their temperature decreases which results in a higher 7. We predict therefore that, in the context of the experiments 
reported in Q, a mixture of steel and aluminum {asteei ~ 0.9, aai « 0.83, nisteei ~ 3TOa/) should yield a smaller value 
of 7 than the brass-glass mixture {ttbrass ~ 0.8, Ugiass ~ 0.9, rubrass ~ 3mgiass) for which the measured 7 is close 
to 0.6 — 0.7. The dependence of 7 upon number fraction Xi = pi/ p may on the other hand be counter-intuitive: at 
a given mean density po, an increase of the relative fraction xi of heavy particles leads to an increase of 7 when the 
heavy particles are the more elastic (see Fig. f6|) . This effect was also clearly observed for the homogeneously heated 
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mixture [|6). On the other hand, an increase of xi leads a relatively weak decrease of 7 when the heavier particles 
are the less elastic, whereas the opposite (albeit also quite weak) trend could be observed in |]l6| . 



7 7a: 7r a 7 7a; Ir 



-1 0.88 0.92 - -1 0.79 0.845 - 

-0.5 0.825 0.89 0.83 -0.5 0.7 0.78 0.69 

0.79 0.86 0.8 0.65 0.74 0.66 

TABLE I: Values of 7, ^y^, 7^ in the middle of the system for A'^ — 500, atj = 0.85, ?7i_o = t?2,o, "^i ~ 3m2 (left) and mi = 5m2 
(right) 



Q 7 7j: 7r Q 7 7a: 7r 

-1 0.735 0.775 - -1 0.95 1. 

-0.5 0.69 0.735 0.735 -0.5 0.89 0.99 0.84 

0.665 0.72 0.72 0.85 0.96 0.81 

TABLE XL Values of 7, 7^, 7^ in the middle of the system for = 500, = 0.9, 0.8, 0.7 (left) and aij = 0.7, 0.8, 0.9 (right), 
mi = 3m2, rjifi = ?)2,o. 



C. Velocity distributions 



As in the monodisperse case, we have measured the single-particle velocity distributions, which are anisotropic 
as expected. The vertical velocity distributions are similar to those shown in Fig. ||, and the horizontal velocity 
distributions show strong non-Gaussian features, as in the monodisperse case. Moreover, it appears in Fig. |l^ that 
the rescaled velocity distributions Pi{cx) and P2{cx) are very close (even if not equal, see also [^) for both types 
of particles. The differences between Pi(cx) and P2{cx) increase if the inelasticities or the mass ratio increase. 
Pi{cx) depend slightly on the various parameters, in the same way as the velocity distributions of the monodisperse 
situation; this dependence would probably be very difficult to measure in an experiment, which would probably lead 
to the conclusion that Pi{cx) ~ P2{cx) ■ 



V. CONCLUSION 



In this study, we have considered vibrated granular gases well outside the Boltzmann limit of (very) low densities. 
The molecular dynamics simulations performed are free of the approximations underlying the usual kinetic theory or 
hydrodynamic approaches. Taking due account of the first correction to the ideal gas contribution in the equation 
of state (second virial order), we however found a remarkable constant yy component of the pressure tensor over the 
whole cell, for monodisperse or bidisperse systems, despite the strong density and temperature heterogeneities due to 
the realistic energy injection mechanism. 

The study of the velocity distributions along the horizontal direction (perpendicular to the energy injection) has 
revealed non-Gaussian features similar to experiments, which depend weakly on the various parameters involved in 
the model. 

The projection of the dynamics onto the horizontal direction has allowed us to gain insight into the correlations 
between the effective restitution coefficient aid and the relative velocities gx of colliding particles. The measured 
conditional probability distributions ^{aid\gx) are in agreement with the forms proposed in based upon partial 
experimental data. The link between ^(aid\gx) and the velocity probability distribution functions JTo[ has been 
confirmed. 

In the case of binary mixtures we have analyzed the ratio of granular temperatures as a function of the various 
parameters, and found a very good qualitative agreement with experiments. The velocity distributions of the two 
components have moreover been shown to be very similar. 
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APPENDIX A: INCLUSION OF A TANGENTIAL RESTITUTION COEFFICIENT 

In this appendix we give the coUision rules when a tangential restitution coefficient is introduced (see also [pl|). 
The two colliding particles, labeled (1) and (2), have masses uii, diameters Ui, moment of inertia = rriiqai/A with 
q = 1/2 for disks and 2/5 for spheres. The precoUiding velocities are Vi,uJi, and postcoUiding velocities are denoted 
with primes. 

The normal unit vector is defined as: 

^ = f^. (Al) 
\ri - rzl 

The relative velocity of the contact point 

g = vi - V2 - ( Y'^i + T'^^) ^ ^ 
has normal component Qn = {o ■ 5^)5^ and tangential component Qt = 9 Qn (this defines the tangential unit vector 

* = gtl\gt\- 

The postcollisional velocities can be expressed simply in terms of the precollisional velocities through the introduc- 
tion of the linear momentum change of particle (1) 

l^P ^ mi{v[ ~ vi) ^ -m2{v2 - V2) . (A3) 
Indeed the change of angular momentum is 

— {oj[-u}i)^-^x AP (A4) 

One obtains: 

AP 

v[ = Vi-\ (A5) 

mi 

AP 

V2 = V2 (A6) 

1712 

ujI = LJ, -^axAP (A7) 

The normal and tangential components of AP are then computed using the definition of the normal and tangential 
coefficients of restitution: 



£/; = (A8) 
g't - • (A9) 



Since g„ = [{vi — V2) ■ a-] a, the first relation leads to 



AP.? = {I + a){vi ~ V2) ■ a . (AlO) 

nil + W2 



Using the definition of gt, and with li — niiqai/A, one obtains also 



[where APt = {AP ■ Finally, 



g'^=g, + APt{— + —)(l + -] (All) 
nil 1712 J \ q 



mim2 ( , 1 + a* 



AP = ^ (1 + + T-T^9* (A12) 

TOi + TO2 \ 1 + 1/g 
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APPENDIX B: EQUATION OF STATE FOR A POLYDISPERSE INELASTIC MIXTURE 



In this appendix, we adopt a kinetic definition of the total pressure and compute this quantity for an arbitrary 
homogeneous mixture of species i, with number fraction Xi = Pi/p. Invoking the virial theorem, the excess pressure 
pcx — p _ pidcai ~ p _ J2i£iTi is related to the collisional transfer of linear momentum: the partial excess pressure 
of species i reads (see e.g. pj) 



E 



lim — — - / I J, 

coll. partner of i 

lim 

t^oo dV t ^ 



j. coll. partner of i 



rrii + m 



■(1 + aij)(d- ■ Vij) Gij where a. 



(Bl) 
(B2) 



In these equations, it is understood that the summation runs over all the collision events involving a particle of type i 
and an arbitrary partner j, in a large volume of measure V . The collisional transfer appearing in Eq. (B2) is readily 



computed within Enskog-Boltzmann kinetic theory, where the velocity distribution functions '^i{v) obey the set of 
non-linear equations 



dt(piivi,t) = ^ XijO-fj ^TLj j dv2 / da-Q{a ■ Vu) {a- ■ Vu) 
j=i J J 



(B3) 



where 9 denotes the Heavyside distribution and (y\^vl2) are the pre-coUisional velocities converted into (vx^v^) by 
the collision rule (|l|)-(0). Equation (B2) may be rewritten 



, A/'a „ „ 

Pr = zZ^^i'^ii^ / ^■"i^^2 / daQ{a ■Vi2){^ ■Vi2)ipi{vi)ip.j{v2) ^— ^(1 + ay)(ff • ^12) cTij. (B4) 

Zu . J J TTli -j- JJlj 



Summing the contributions of all species, the total excess pressure follows: 



pc 



2d 



2d ^ 



ij n^ rij 



mtirij ^-^^ I dvidv2 I da Q {a ■ V12) {a ■ Vi2)'^ (fi{vi)(pJv2) (B5) 
rrii + rrij ' ' 



rrii nij 
rrii + Tij 



(l + ay) 



J daQ {a- ■ V12) {a- ■ v^f J dvidv2 {v^ + V2) ifi{vi)ipj{v2){B6) 



where V12 is the unit vector along V12 , and where the contribution from the dot product Vi ■ V2 vanishes by symmetry 
in the last integral. The integral over the solid angle is related to the volume Vd of a sphere with diameter 1: 



j daQ{a ■ V12) (ct • V12 



rd/2 



dT{d/2) 



(B7) 



where T is the Euler function and it is understood that V12 denotes an arbitrary unit vector in (B7). The volume 
Vd is itself related to the packing fraction 77 through 77 — pVd{cr'^). From the definition of kinetic temperatures 
/ v'^(pi{v) dv = dTi/rrii, we get 



mi m„ 



(B8) 



from which we deduce the equation of state (|4|). In this last step, no approximation (e.g. Gaussian etc) is made 
concerning the Lpi. On the other hand, the computation of any other moment [a ■ V12Y than p = 2 requires the 
det aile d knowledge of the velocity distributions ||2^ . It is also noteworthy that the decoupling of velocities Vi and V2 
in (B6) is a specific property of the momentum transfer, which significantly simplifies the calculation. 



APPENDIX C: HYDRODYNAMICS 



In this appendix, we recall the hydrodynamical approach considered by Brey et al. 
of energy injection at both boundaries y = and y = L. The situation investigated in 



, and adapt it to the case 
is that of a vibrating wall 
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at 2/ = 0, and a reflecting wall at j/ = L so that the temperature and density gradients vanish aX y = L. In our 
no-flow configuration with two vibrating walls, the gradients vanish by symmetry in the middle of the cell (y = L/2), 
so that restricting to y £ [0, L/2] allows to use directly the expressions derived in |^ (which amounts to the formal 
identification y — > 2y and N — > N/2. For completeness and clarity, we will however adapt the argument to our 
geometry. 

In the case of a stationary system, without macroscopic velocity flow, the hydrodynamic equations reduce to 

V-P = (CI) 

—y-q + TC = . (C2) 
pa 

Here P is the pressure tensor, q is the heat flux, and ( the cooling rate due to the collisional energy dissipation. In 
the Navier-Stokes approximation for a low density gas described by the Boltzmann equation modified to account for 
the inelastic nature of collisions 

P = PI (C3) 
q = ~K\7T-fi\7p (C4) 

where P is the ideal gas pressure: P = pT. The explicit expressions of the heat conductivity k, the transport 
coefficient p, and cooling rate C may be found in j^]. The important ingredient is that p is proportional to T'^^'^ / p and 
K to VT, while C cx p/\/T^ with coefficients depending on the inelasticity a. 

The system is considered homogeneous in the x direction, so that only gradients along the y direction are taken 
into account. We emphasize that the ideal gas equation of state P = pT is assumed, and this simplification is an 
important ingredient in the following derivation. The previous equations then reduce to: 

dP 

m.)3 (,^m-\_^^ (C6) 



dp dy \ dy 

In order to simplify the equation on the temperature, it is convenient to introduce a new variable ^ defined by 

dy 



d^ - vM^tA = Ca''-'y/a{a)p{y)dy (C7) 

Ky) 

where A(j/) = [Ca"^^^ p{y)]^^ is the mean-free-path (C — 2\/2 for d ~ 2), and a{a) includes all the dependence in a. 
Equation (C6) now reads 

^Vf=Vf. (C8) 

The variable ^ takes values between and ^m, with cx N. Then ^/T = Aexp(— ^) -I- i3exp(^) where A and B 
depend on the boundary conditions. In the case of two vibrating walls, the solution is symmetric with respect to 
y = L/2 (or ^ = Cm/2). With T(0) = T(C^) = Tq one obtains 

TiO = . ll, (sinh(U - + sinh if . (C9) 
smh i„i 



It is possible to integrate d^ = Ca'^ ^ \J a{a)n{y)dy — Ca*^ ^ a{a)pdy /T [y) to obtain y{i) and P: 

y_ ^ C + sinhCcosh(C^ ~0 
L i,n + smhi,n 

Those equations are the same as for the case of one vibrating wall Q, but with — > 2^™ and i — > 2i, as expected 
on the basis on the symmetry argument proposed above. It is possible to invert T{i) and therefore to obtain the 



11 



profiles y{T) (two symmetric branches): 



^ ± cosli ^ 



^ 1 

— cosli — 
Jo ^ 



y _ C + sinli^cosh($„t - Q 
L + sinh 



(C12) 
(C13) 
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FIG. 1: Density profiles for two normal inelasticities and two densities. In all cases, the number of particles is A'' = 500. The 
symbols correspond to the smallest density (the mean packing fraction, averaged over the whole system is tjo = 0.015) and the 
lines are for a higher density (770 = 0.04). The ratio r?(2/)/r?o is also the ratio p{y)/ po of local density normalized by the mean 
one. 
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FIG. 2: Density profiles r){y)/r]o (upward curves) and temperature profiles (downward curves) for a given normal restitution 
coefficient a = 0.9, and different tangential restitutions (AT = 500 particles, mean packing fraction 770 = 0.015). The temperature 
is the total one (including horizontal and vertical degrees of freedom); it is expressed in arbitrary units but all curves correspond 

to the same velocity of the vibrating piston. From top to bottom for the temperature T{y) and from bottom to top for the 
density, the curves correspond respectively to a* = — 1, a* = —0.8, a* = —0.5 and a* = 0.2. 
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FIG. 3: Temperature profile for a = 0.9 and ?7o = 4%. Tlie liorizontal T^, vertical Ty and total temperature T — {T^ +Ty)/2 
are shown. Inset: anisotropy factor A = {Ty — Tx)/{2T) as a function of height. 
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FIG. 4: Pressure given by the equation of state (^). (a) The symbols correspond to P = p(y)T{y)[l + (1 + ct)'q{y)x{y)] (see 
text), where T is the total temperature. The lines immediately below a given set of symbols show the ideal gas contribution 
p{y)T{y) only. For the three situations investigated, the mean density is the same (t^o = 0.04). 

(b) Same figure with the vertical temperature Ty instead of T inserted in the equation of state, yielding therefore the yy 
component of the pressure tensor. 
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y/L 

FIG. 5: Fits of the temperature profiles measured in MD witli the analytical expression (Q. The fits are shown with continuous 
curves while the symbols stand for the MD measures. For clarity the fits are restricted to heights L/2 < y < 0.8L. 




FIG. 6: Probability distribution function (p.d.f.) of the vertical velocity component Cy — Vy/y^Ty, for different heights. By 
definition, (Cy) — 1 whatever the altitude y. Here, rjo — 0.04, A'^ = 500, a = 0.9 and a* = 0. 
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FIG. 7: (a): Probability distribution function of the rescaled horizontal velocity component Cx — Vx/VT^, on a linear-log plot. 
Here rjo = 0.015, A'^ — 500, a = 0.9 (pluses) and 0.8 (stars), and a' = 0. The solid line is the Gaussian with unit variance, the 
circles correspond to experimental data [0, ^ for steel beads. 

(b): Probability distribution function of the angular velocities for the same parameters. 
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FIG. 8: Probability distribution function of effective one-dimensional restitution coefficients a^d- The MD results are compared 
to the experimental measures of Feitosa and Menon Q on steel and glass samples (for which the nominal restitution coefficient 
may be considered close to 0.9). 
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FIG. 9: (a): Conditional p.d.f. of aid for a given value Qx of order unity. Note the different shapes for thermal and vibrating 
walls. 

(b): Same, as a function of {aidQx)^ (and Qx = 0.2, 0.5, 1., 1.5, 2., 3., 4., 5.) for vibrated walls with a = 0.9, q* = and t^o = 0.015. 




FIG. 10: Probability distribution function of energy restitution coefficients f3. Various tangential restitution coefficients a* are 
considered for a — 0.9 and 770 ~ 1.5%. The circles represent the experimental data for steel grains [w] 
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FIG. 11: (a): The symbols show the pressure calculated from the complete equation of state for a binary mixture ( |ll|) including 
Enskog correction, while the lines immediately below display the ideal gas contribution pi(y)Ti(y) + p2{y)T2{y) to the pressure. 
The three sets of curves correspond to: upper set rjo = 0.015, an = 0.9, ai2 ~ 0.8, Q22 = 0.7, mi — 5m2; middle set rjo = 0.04, 
ail — 0.9, ai2 = 0.8, 022 = 0.7, mi = 3m2; lower set r)o = 0.04, an — 0.7, ai2 — 0.8, a22 ~ 0.9, mi = 3m2. 
(b): same curves, where the temperatures are the vertical ones Ti,y instead of the total Ti — {Ti,x + Ti,y)/2, yielding therefore 
the yy component of the pressure tensor. 




FIG. 12: Temperature profiles for an equimolar granular mixture, driven by vibrating walls. The symbols show the MD 
measures, and the lines are fits to the analytical expression derived for the single component case. In all cases, the particle 1 
(the heaviest) has mass mi — 3m2; its temperature Ti corresponds to the two lower sets. 
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FIG. 13: (a): Vertical profiles for a binary mixture with mi = 3m2, r?o — 0.015, and equal mean densities 771,0 ~ ?72,o (excitation 
by vibrating walls). From bottom to top: temperature profiles of both species, density profiles 7/2(2/) /(2r;o) 771(7/)/ (2770). 
Since ai = (T2, the packing fraction rii is proportional to the local density pi of species i. The upper dashed curve shows the 
temperature ratio 7 = T2/T1 as a function of height, and the circles show the same quantity for a non equimolar mixture where 
Vi,o = 8772,0- 

(b): Same with a higher mass ratio 77ii — 67712. 
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FIG. 14: Density profiles and temperature ratio profiles (binary mixture, vibrating walls). The lines correspond to aij = 
0.7; 0.8; 0.9 whereas the symbols are associated with "reverse" inelasticities aij — 0.9; 0.8; 0.7. The other parameters are a* = 0, 
7711 = 37772, 7/1,0 = f?2,o, rjo = 27/1,0 ~ 0.015. The upper fiatter curves (dashed line and stars) show the temperature ratio. As in 
Fig. |l^, the density of heavy particles pi (thick continuous curve and circles) is more peaked and denser in the middle of the 
cell than that of light grains (thin continuous curve and squares). 
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FIG. 15: Effect of density on the temperature ratio for mi = 3m2, ciij = 0.9; 0.8; 0.7 (vibrating walls). Graph (a) shows 
the total ratio T-z/T-i and graph (b) shows the ratio of horizontal temperatures T2,x/T\^x- In both cases, the corresponding 
experimental measures are shown in the insets for a steel glass mixture (at different densities, but with a density ratio of 2, 
close to that of the MD simulations 0.04/0.015 ~ 2.6). The purpose is to show that the changes induced by density in MD are 
qualitatively the same as in the experiments. 
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FIG. 16: Influence of number fraction on the temperature ratio T2/T1. The total number of particles is N = Ni -\- N2 = 500 
(vibrating walls). Given that ai = (72, N1/N2 = 8 corresponds to 7ji,o = 8772,0- 
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FIG. 17: Probability distribution functions of the rcscalod horizontal velocity components ct^x = Vi,x/ y/TXx, for an cquimolar 
mixture. Squares are for Pi (heavy grains) and circles for P2 (light grains). Here rjo = 0.015, A'' = 500, aij = 0.9,0.8,0.7, 
mi = 3m2 and a* = 0. The solid line is the Gaussian with variance 1. 



